#include "finite_diff_method.h"
#include "head.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h> // Get_Time
#include <sys/resource.h> // Get_Time
#include <mpi.h>


/**
 * @file   finite_diff_method.cpp
 * @author W. P. Yao <ywp@pku.edu.cn>
 * @date   2015年11月12日15:54:22
 * brief   有限差分方法
 */

extern int rank, nprocs;

int N;                    // 网格剖分参数
int Points_Count;         // 网格节点总数
 
int *Points_Class;        // 按进程为节点分类
int *Points_Class_Count;  // 每组节点数目

int *Points_Class_Min;    // 计算时用到点的最小指标
int *Points_Class_Max;    // 计算时用到点的最大指标
int Min_P, Max_P;         // 数据散发

int *Points_Coordinate;   // 节点坐标信息，三列表
int *Points_Zone;		  // 节点所属区域或边界的信息，一维数组
int *Table_Connection;    // 节点连接信息，六列表

double Mesh_Size;		  // 空间步长
double Time_Step_Ratio;   // 时间步长比
double Time_Step;         // 时间步长
double Ratio;             // 网格比 tau/h^2
int Step_Max; 			  // 最大计算步数

/**
 * 从文本文档中读取网格数据
 * @return 读取成功返回0，否则返回1/2
 */

int Read_Mesh()
{
	FILE *fp1 = NULL;
	FILE *fp2 = NULL;

	fp1 = fopen("../Data/Point_Table.txt","r+");
	fp2 = fopen("../Data/Connection_Table.txt","r+");

	if(fp1 == NULL || fp2 == NULL) return 1;

	FILE *fp0 = NULL;
	
	fp0 = fopen("../Data/Mesh_Info.txt","r+");
	
	if(fp0 == NULL) return 2;

	fscanf(fp0,"%d\n%d\n", &N, &Points_Count);

	fclose(fp0);
	fp0 = NULL;

	Mesh_Size = 1.0/N;     // 空间步长

	Points_Coordinate = (int*)malloc(Points_Count*3*sizeof(int));
						  // 坐标
	Points_Zone       = (int*)malloc(Points_Count  *sizeof(int));
						  // 区域或边界标记
	Table_Connection  = (int*)malloc(Points_Count*6*sizeof(int));
						  // 节点连接表   ???

	
	int AN;
	for(int m=0;m<Points_Count;m++){
		fscanf(fp1,"%d %d %d %d %d", &AN,
			Points_Coordinate+3*m,
			Points_Coordinate+3*m+1,
			Points_Coordinate+3*m+2,
			Points_Zone+m);

		fscanf(fp2,"%d %d %d %d %d %d %d", &AN,
			Table_Connection+6*m,
			Table_Connection+6*m+1,
			Table_Connection+6*m+2,
			Table_Connection+6*m+3,
			Table_Connection+6*m+4,
			Table_Connection+6*m+5);
		}
	
	printf("\n\ttest: %d\n", Table_Connection[6*AN]);
						   // 测试读取是否正确
	printf("Enter to continue...");
	getchar();


	fclose(fp1);
	fclose(fp2);
	fp1 = NULL;
	fp2 = NULL;
	
	return 0;
}

/**
 *
 * 并行的读取网格数据
 * 由进程0读和广播
 * @return 读取成功返回0，否则返回1
 */

int MPI_Read_Mesh() {
	
	if (0 == rank) Read_Mesh();

	MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&Points_Count, 1, MPI_INT, 0, MPI_COMM_WORLD);

	Mesh_Size = 1.0/N;

	if (0 != rank){
		
		Points_Coordinate = (int*)malloc(Points_Count*3*sizeof(int));
		Points_Zone = (int*)malloc(Points_Count*sizeof(int));
		Table_Connection = (int*)malloc(Points_Count*6*sizeof(int));

		}
	
	MPI_Bcast(Points_Coordinate, Points_Count*3, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(Points_Zone, Points_Count, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(Table_Connection, Points_Count*6, MPI_INT, 0, MPI_COMM_WORLD);
	
	// 节点分组
	Points_Class = (int*)malloc((nprocs+1)*sizeof(int)); // 节点分组指标
	Points_Class_Count = (int*)malloc(nprocs*sizeof(int));
	
	int cn = Points_Count/(nprocs-1);      // 商
	int rn = Points_Count - cn*(nprocs-1); // 余数
	
	// 进程0不做计算
	Points_Class[0] = 0;
	Points_Class[1] = 0;
	Points_Class_Count[0] = 0; 

	for(int i=1; i<nprocs; i++){
		
		if(i<=rn){
			Points_Class[i+1] = Points_Class[i] + cn + 1; // 连续分组，每组个数尽量相同
			Points_Class_Count[i] = cn + 1; // 每组节点个数
			}
		else{
			Points_Class[i+1] = Points_Class[i] + cn;
			Points_Class_Count[i] = cn;
			}
		}
	
	Max_P = 0;
	Min_P = Points_Count;
	for(int i = Points_Class[rank]; i < Points_Class[rank+1]; i++){
		for(int j = 0; j < 6; j++){
			int P = Table_Connection[6*i+j];
			if(P < 0) continue;
			Min_P = Min_P>P?P:Min_P;
			Max_P = Max_P>P?Max_P:P;
			}
		} // 每组依赖的数据长度
	
	if(0 == rank){
		Min_P = 0;
		Max_P = -1;
		}
	Points_Class_Min = (int*)malloc(nprocs*sizeof(int));
	Points_Class_Max = (int*)malloc(nprocs*sizeof(int));

	MPI_Allgather(&Min_P, 1, MPI_INT, Points_Class_Min, 1, MPI_INT, MPI_COMM_WORLD);
	MPI_Allgather(&Max_P, 1, MPI_INT, Points_Class_Max, 1, MPI_INT, MPI_COMM_WORLD);

	printf("rank=%d:\t%d~%d\nused:\t%d~%d\n", rank, Points_Class[rank],
		Points_Class[rank+1], Min_P, Max_P);

	return 0;
	}


/**
 * 设定初值
 * @param u 设定初值
 */

void Init_Val(double *u){
	
	for(int i=0;i<Points_Count;i++){
		*(u+i) = 0.5;  
		}
	}

/**
 * 设置计算时间步长
 * @param t0, 初始时刻
 * @param r,  时间步长比
 * @param n,  最大计算步数
 */

void Set_Time_Step(double t0, double r, int n){
	
	Time_Step_Ratio = r;
	Step_Max = n;
	}


/**
 * 自定义数据散发函数
 * 进程之间相互补充需要的数据
 *
 * @param u 要散发的数据
 *
 * @return 0
 */
int MPI_My_Scatter(double *u){
	
	MPI_Status status;
	int left  = (nprocs+rank-1)%nprocs;
	int right = (rank+1)%nprocs;

	MPI_Sendrecv(u+Points_Class[rank],
			     Points_Class_Count[rank],
				 MPI_DOUBLE, left, 99,
				 u+Points_Class[right],
				 Points_Class_Count[right],
				 MPI_DOUBLE, right, 99,
				 MPI_COMM_WORLD, &status); // 向左边传数据
		
	MPI_Sendrecv(u+Points_Class[rank],
			     Points_Class_Count[rank],
				 MPI_DOUBLE, right, 100,
				 u+Points_Class[left],
				 Points_Class_Count[left],
				 MPI_DOUBLE, left, 100,
				 MPI_COMM_WORLD, &status); // 向右边传数据
	
	//printf("\t rank=%d, send over\n", rank);
	
	return 0;	
	
	}



/**
 * 逐步计算
 */

void Run(){
	
	double tau = pow(Mesh_Size,2)/6;   // ??? 
	Time_Step = Time_Step_Ratio * tau; // ??? 
	
	if(0 == rank){
		printf("\tThere are %d points for N = %d.\n", Points_Count, N);
		printf("\tMesh_Size = %f, Time_Step_Ratio = %f\n", 
			Mesh_Size, Time_Step_Ratio);
		printf("\ttau = %f\n", tau);
		printf("\tTime_Step = %f\n", Time_Step);

		printf("Enter to start...");
		getchar();
		}

	double *u1 = (double*)malloc(Points_Count*sizeof(double));  // 数值解i
	double *u2 = (double*)malloc(Points_Count*sizeof(double));  // 数值解i+1
	double *u3 = (double*)malloc(Points_Count*sizeof(double));  // 数值解vs

	Init_Val(u1);
	Init_Val(u3);
	
	if(0 == rank) Write_Data(0,u1);   							// 写数据@t=0

	double cpu0, cpu1, wall0, wall1;
	Get_Time(&cpu0, &wall0);

	double *u_old = u1;
	double *u_new = u2;

	Ratio = Time_Step/pow(Mesh_Size,2); 						// 网格比???

	int Step = 0;
	while(Step++ < Step_Max){
		//for(int m=0;m<Points_Count;m++){                      // 串行
		for(int m=Points_Class[rank]; m<Points_Class[rank+1]; m++){	
			switch(Points_Zone[m]){
				case NORM:
				Normal_Point(u_old, u_new, NORM, m); break; //正常节点
				case INNER:
				Inner_Point(u_old, u_new, INNER, m); break; //内边界节点
				case OUTER:
				Outer_Point(u_old, u_new, OUTER, m); break; //外边界节点
				default:
				break;
				}
			}

		if(Step%1000 == 0){                                  //隔1000步输出一次
			
			MPI_Gatherv(u_new+Points_Class[rank], Points_Class[rank+1]
						-Points_Class[rank], MPI_DOUBLE, u_new, 
						Points_Class_Count, Points_Class, MPI_DOUBLE,
						0, MPI_COMM_WORLD);
			
			double Norm_Inf = 0;

			if(0 == rank) {
				Write_Data(Step, u_new);

				for(int m=0;m<Points_Count;m++){                 //计算无穷范数
					double d = fabs(*(u3+m) - *(u_new+m));
					if(d > Norm_Inf)
						Norm_Inf = d;
						*(u3+m) = *(u_new+m);
					}
				printf("Step %d:\tdiff = %5e", Step, Norm_Inf);
				
				Get_Time(&cpu1,&wall1);

				printf(",\tTakes %f second\n",wall1-wall0);
		    	}
			
			MPI_Bcast(&Norm_Inf, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
		
			//if(Norm_Inf < 1e-9) break;   // 终止条件
			}
		
		MPI_My_Scatter(u_new);

		double *temp;    // 交换存储空间
		temp  = u_new;
		u_new = u_old;
		u_old = temp;
		}
	
	MPI_Barrier(MPI_COMM_WORLD);
	free(u1);
	free(u2);
	
	}



/**
 * 正常节点的计算，6个相连点都是有效地
 *
 * @param u 上一时刻的温度
 * @param u_new 此时的温度
 * @param flag  区域的标记
 * @param m     节点的编号
 */

void Normal_Point(double *u, double *u_new, int flag, int m){
	
	int locate = 6*m;
	double u_OOO = *(u+m);
	double u_MOO = *(u + Table_Connection[locate++]);
	double u_POO = *(u + Table_Connection[locate++]);
	double u_OMO = *(u + Table_Connection[locate++]);
	double u_OPO = *(u + Table_Connection[locate++]);
	double u_OOM = *(u + Table_Connection[locate++]);
	double u_OOP = *(u + Table_Connection[locate++]);

	double LAPLACE_u = u_MOO + u_POO + u_OMO + u_OPO + u_OOM
		+ u_OOP - 6*u_OOO;
	*(u_new+m) = u_OOO + Ratio*LAPLACE_u;
	} 

void Inner_Point(double *u, double *u_new, int flag, int m){
	*(u_new+m) = 0.0;
	}

void Outer_Point(double *u, double *u_new, int flag, int m){
	*(u_new+m) = 1.0;
	}


/**
 * 将第n步的计算结果写入文本文件
 * 
 * @param n 步数
 * @param u 温度
 * 
 * @return 写入成功返回0，否则返回1
 */

int Write_Data(int n, double *u){
	
	FILE *fp = NULL;
	char Name[100];
	sprintf(Name,"../Result/Step %d.txt", n);
	fp = fopen(Name,"w+");
	if(fp==NULL) return 1;

	for(int m=0;m<Points_Count;m++)
		fprintf(fp,"%f\n",*(u+m));
	
	fclose(fp);
	fp = NULL;

	return 0;
	}


/**
 * 计算程序运行的时间
 *
 * @param cpu  CPU时间
 * @param wall 墙上时间
 */

void Get_Time(double *cpu, double *wall){
	
	struct timeval tv;
	struct rusage  ru;               //resource.h
	if(cpu!=NULL){
		getrusage(RUSAGE_SELF,&ru);  //获取当前进程的资源使用信息
		*cpu = ru.ru_utime.tv_sec +
			(double)ru.ru_utime.tv_usec*1e-6; //user time used
		}
	if(wall!=NULL){
		gettimeofday(&tv,NULL);
		*wall = tv.tv_sec + (double)tv.tv_usec*1e-6;
		}

	} 

/**
 * 释放内存
 */


void Finalize(){
	
	free(Points_Coordinate);
	free(Points_Zone);
	free(Table_Connection);
	
	free(Points_Class);
	free(Points_Class_Count);

	free(Points_Class_Min);
	free(Points_Class_Max);
	}



